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■ ■ We expand the Potts-percolation model of a solid to include stress and strain. 

Neighboring atoms are connected by bonds. We set the energy of a bond to be given 

j/3 \ by the Lennard-Jones potential. If the energy is larger than a threshold the bond is 

C^ ■ more likely to fail, while if the energy is lower than the threshold the bond is more 

likely to be alive. In two dimensions we compute the equation of state: stress as 

C \ function of inter-atomic distance and temperature by using renormalization group 

^ ' and Monte Carlo simulations. The phase diagram, the equation of state, and the 

CSJ ■ isothermal modulus are determined. When the Potts heat capacity is divergent the 

t^^ , continuous transition is replaced by a weak first-order transition through the van der 

m 

Waals loop mechanism. When the Potts transition is first order the stress exhibits 
t~^ ' a large discontinuity as function of the inter-atomic distance. 

o 

PACS numbers: 05.10.Ln,05.10.Cc,62.20.-x 

I. INTRODUCTION 

The mechanical properties of solids, such as melting^ '^ and mechanical failure^"-, are 
topics of considerable interest. In this paper we continue the analysis of an equilibrium 
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statistical mechanics inodeF'^ of a solid. Previously we assumed^ harmonic springs and 
evaluated the role of thermal fluctuations by using renormalization-group and Monte-Carlo 
simulations. Furthermore we studied^- the model with an extended defect line and found 
a hybrid, first- and second-order phase transition. In this paper we study the equation 
of state of the solid, stress as function of strain and temperature, obtained by assuming 
neighboring atoms are separated by a fixed inter-atomic distance and associating to each 
pair of neighboring atoms the Lennard- Jones energy. An alternative realistic anharmonic 
energy versus atomic distance due to Ferrante^i has been considered in a previous work—. 

The model is defined in Section II. We assume the energy of a pair of neighboring atoms 
to be given by the Lennard- Jones 6-12 potential. If the energy of such a spring is larger 
than the threshold energy, the probability for its failure is higher than 50%. This model is 
mapped into a Potts model with couplings that are dependent on the inter-atomic distance. 
The free energy, number of live bonds, their fluctuations, stress, and modulus are computed 
using renormalization group and Monte Carlo techniques. Fixing the inter-atomic distance 
to a value independent of the atoms locations makes this model mean-field like. The role 
of the inter-atomic distance fluctuations will be studied in a future work. The number of 
Potts states q is a fugacity conjugated to the number of clusters of live bonds, with g = 1 
corresponding to springs failing independently of each other. A challenging question not 
addressed here is how to connect q to data from a real solid. 

The equation of state is studied paying particular attention to the solid failure signaled 
by an extremum in the dependence of stress on the inter-atomic distance. Beyond the 
maximum under expansion and the minimum under compression for stress versus strain 
dependence, the solid ceases to be thermodynamically stable. The phase diagram includes 
the Potts transition line. If the Potts transition line is in the stable region of the phase 
diagram and if the q value is such that the Potts heat capacity is divergent, a remarkable 
phenomenon occurs in the vicinity of the transition line. A van der Waals loop^ developes in 
the stress-strain dependence signaling a weak phase transition that replaces the continuous 
Potts phase transition. While in the renormalization group calculationsi^ii^ that are exact^ 
on hierarchical lattices^^'^'^the Potts transitions are always continuous, in the Monte Carlo 
simulations we can see both continuous transitions (for small q) and discontinuous transitions 
(for large q). This allows us to explore the influence of the order of the Potts transition on the 
mechanical properties of the solid, such as stress dependence on temperature and interatomic 



distance. 

In Section III we present numerical results based on the renormalization-group Migdal- 
Kadanoff scheme. Monte-Carlo simulations are presented in Section IV. Our concluding 
remarks are found in Section V. 

II. MODEL 

The energy of any pair of neighboring atoms is: 

E{r) = e[Cjr-2Cff] (1) 

where r is the inter-atomic distance and tq is the equilibrium inter-atomic distance (under 
zero stress). If the energy of the spring is larger than the threshold energy Eq the bond is 
more likely to fail than to be alive, p is the probability that the bond is alive and 1 — p the 
probability that the bond is broken. We assume the probabilistic weight w = p/{l — p) io 
depend on energy through the Boltzmann weight: 

w = — ^— = e >=bt (2) 

1-p ^ ^ 

We allow for correlations between failing events by using the Potts number of states q, 
which plays the role of a fugacity controlling the number of clusters. For the same number 
of live bonds, graphs with more clusters are favored if g > 1, while if g << 1 there is a 
tendency to form a few large clusters. If g = 1 bonds fail independently of one another, i.e. 
random percolation process. 

The partition function is obtained^^ by summing over all possible configurations of bonds 
arranged on the lattice 

z= y: f^^ (3) 

C is the number of clusters, including single site clusters, and B is number of live bonds. 

The free energy per bond is / = In Z/Nbonds- The derivatives of the free energy / with 
respect to w and q provide respectively the number of live bonds b and the number of clusters 
c, each normalized by the total number of lattice bonds 



' = % (') 

The derivatives of b and c with respect to w and q provide the fluctuations (variances) 
of those quantities: 

Ab^ = w^— (6) 

ow 

The stress a is calculated by taking the derivative of the free energy with respect to r, 
the inter-atomic distance: 

a = -ksT^ = b—^ 8 

or ar 

The equation of state, Eq. ([H]), states that the average stress is equal to the stress 
associated with each live bond dE{r)/dr multiplied by the number of live bonds b. The 
energy gradient is obtained from Eq. ([1]): 



dE{r) 12e 



:7)"-(^)' 



(9) 



dr tq 

The equation of state yields the isothermal linear modulus (inverse compressibility) m: 

(11) 

We perform Monte Carlo simulations of the model using its mapping into the Potts 
model. By using the Kasteleyn-Fortuin expansion22i2i we can rewrite the partition function, 
Eq. (ED, as: 



The Hamiltonian is: 



Tr^e "bt (12) 



^ Y: J{r)6{s,,s,) (13) 
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where s, is a Potts^^ spin at the lattice site i taking q values. The coupling constant J{r) is 
related to the original parameters by 



E{r)-Eo 

J{r) = ln(l + w) = ln(l + e "bt ) (14) 

To get the equation of state from Monte Carlo simulations for a given T, r, and Eo, we 
calculate the Potts coupling constant J(r), using Eq. f lT^ . Then we rewrite Eq. (|8]) to get 
the stress. The number of live bonds b is obtained from the Potts energy 



u = -6{si,Sj) (15) 



The equation of state becomes 



, Eir)-E,, ,dE(r) , , 

a = -u\e ''bt + 1-1^^ (16) 

dr 

In all numerical results that will be presented next in Section III from renormalization 
group calculations and in Section IV from Monte Carlo simulations, we express energy in 
units of e, temperature in units of e/A;^, distance in units of vq and stress in units of e/ro- 

III. RENORMALIZATION GROUP 

The Migdal-Kadanoff— »^ recursion equation for two dimensions is 

w'=[l + -^]^ - 1 (17) 

The free energy / = In Z/Nb, Nb being the number of lattice edges, is 

/ = E^ (18) 

where 

C = 2\n{2w + q) (19) 

The recursion equations (11 71) - (11 9 j) represent the exact solutions-^ for the diamond hierar- 
chical lattice^ii^. 

The renormalization group flows are governed by the following flxed points: i. w = 
(non-percolating live bonds), ii. w = oo (percolating network of live bonds), iii. w = Wc 
(Potts critical point). Using the free energy we can compute the number of live bonds b, their 
fluctuation (variance), stress, thermal expansion and modulus. Each of those quantities is 
scaled by the total number of lattice bonds. 
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FIG. 1: Phase diagram (T,r) plane for Eq = —0.5, q = 10. Line with crosses is instabihty line 
and line with diamonds is Potts transition line. The model solid is thermodynamically unstable 
outside the instability line. 



For given values of Eq and q, the phase diagram in the plane (T, r) includes two type 
of singularity lines. In Fig. 1 we fixed the threshold energy at Eq = —0.5 and the clusters 
fugacity at g = 10. The Potts- percolation line separates the region where the probability 
for formation of an infinite cluster of live bonds is non-zero from the region where the 
probability is zero. On this line the probability weight defined in Eq. ([2]) is equal to Wc- On 
the instability line the solid becomes soft, as its modulus vanishes da/dr = 0. The minimum 
of stress versus inter- atomic distance under compression (r < 1) and its maximum under 
expansion (r > 1) provide the limits of stability for the solid. The solid becomes soft at 
these points. Beyond this line the solid is thermodynamically unstable since da/dr < 0. 
Note that a large portion of the Potts transition line (solid, red) is situated in the unstable 
region. 

We show in Fig. 2 the stress versus inter-atomic distance for two isotherms T = 0.2 and 
T = 0.26 respectively. Fig.3 contains zooms close to the Potts transitions that show the 
van der Waals loops. The modulus is shown in Fig. 4 for the same values of the model 



parameters. 

The T = 0.2 isotherm is in the region of phase diagram where the Potts transition is in 
the instabihty region while the T = 0.26 one is in the region where the Potts transition is in 
the stable region. A zoom view on the isotherm T = 0.26 in the neighborhoods of the two 
Potts transitions (see Fig. 3) reveals van der Waals loops-*^^ that yield weak discontinuous 
transitions (small discontinuity) by means of the Maxwell construction. Hence discontinuous 
transitions replace the continuous transition of the Potts model for the hierarchical lattice. 
The modulus and the thermal expansion (Fig. 4) show anomalies at the Potts transition 
that are connected to the divergence of A6^ (or of the heat capacity) in the Potts model 
on the diamond hierarchical lattice for q = 10. The modulus becomes negative, and thus 
the solid is thermodynamically unstable, when A6^ is large enough, see Eq. (8). For the 
square lattice we expect the same phenomenon for g = 2,3,4 where the heat capacity (and 
thus A6^) is infinite at the Potts critical point. This will be confirmed through Monte Carlo 
simulations in Sec IV below. 

In Fig. 5 we show the isotherm where the temperature equals the critical temperature for 
r = 1 (equilibrium inter-atomic distance for zero stress). We also show the corresponding 
modulus. The anomaly apparent in Figs. 2-4 is not present at r = 1 because the derivative 
of E with respect to r [see Eq. (9)] is zero at r = 1 and thus the negative contribution on 
the right hand side of Eq. (10) vanishes. 

To illustrate the dependence of the phase diagram on the parameters q and Eq we show, 
beside the phase diagram for Eq = —0.5, g = 10 of Fig. 1, the phase diagram for Eq = 0.5, 
g = 10 in Fig. 6, and the phase diagrams for Eq = —0.5, g = 1, and for Eq = 0.5, g = 1 in 
Fig. 7. The Potts line in the (T, r) plane is obtained by substituting in Eq. (2) the critical 
value Wc{q), Tq = 1 and using Eq. (1): 

(^Y^-2{^f = EQ-THwM) (20) 

For the self-dual square lattice the critical value-- is Wc{q) = ^/q and this will be used 
in Sec. IV below to verify the accuracy of the Monte Carlo simulations. The instability line 
originates at zero temperature at the two values of r for which E{r) = Eq provided the the 
threshold energy is — 1 < E'o < —0.787. The energy value of -0.787 is the energy at the 
inflexion point r = 1.11. If the threshold energy Eq > —0.787 the instability line in the 
expansion region r > 1 originates from the inflexion point of E{r), while under compression 
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FIG. 2: Equation of state isotherms: stress (white circles) versus strain for Eq = —0.5, q = 10, 
T = 0.2 (upper) and 0.26 (lower). The Potts critical points are indicated by the lines with crosses 
and the stability limits are indicated by the solid vertical lines. The region between the solid lines 
is the stable region. The compression region to its left and the expansion region to its right are 
unstable.. In Fig. 3 we show zooms close to the Potts critical points. 
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FIG. 3: Close-ups on isotherm £^0 = —0.5, q = 10, T = 0.26 near the Potts transition hne under 
compression and expansion, respectively, van der Waals loops, that signal discontinuous transitions 
through the Maxwell construction, are exhibited at the Potts critical points (indicated by vertical 
lines with crosses). The solid line shows the stability limit. 
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FIG. 4: Isothermal modulus m versus inter-atomic distance r for Eq = —0.5, (7 = 10, T = 0.2 
(upper) and 0.26 (lower). The Potts transitions are indicated by the vertical lines with crosses, 
while the solid lines show the stability limits. The region between the solid lines is the stable 
region. The compression region to its left and the expansion region to its right are unstable.. 
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FIG. 5: Stress and modulus versus inter-atomic distance at the critical temperature corresponding 
to r = 1. The Potts critical point is indicated by the line with crosses. No anomaly occurs because 
of the vanishing of -^. The region between the solid lines is the stable region. The compression 
region to its left and the expansion region to its right are unstable.. 
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FIG. 6: Phase diagram (T, r) plane for Eq = 0.5, q = 10. Lines displayed by crosses are instability 
lines and the line with black circles is the Potts transition line. The model solid is thermodynam- 
ically unstable outside the instability lines. 

the starting point is at the r for which E{r) = Eq. 

IV. MONTE CARLO SIMULATIONS 
A. General remarks 



We perform here Monte Carlo simulations of the Potts model governed by the Hamiltonian 
in Eq. (TT3|) with the nearest-neighbor interaction J{r) given by Eq. (TT^ . For a given q, the 
main parameters are T, r and Eq. As in the previous section, we fix q and Eq and make 
simulations in the space (T,r). 

The details of the simulations are the following. We consider a square lattice of size 
N^ X Ny where N^ = Ny = 40, 60, 80, 100. Each lattice site is occupied by a g— state Potts 
spin. We use periodic boundary conditions. Depending on the location of the studied point 
in the phase space, we used an equilibrating time from 10^ to 10^ MC steps per spin and an 
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FIG. 7: Phase diagram {T,r) plane for Eq = —0.5 (upper) and Eq = 0.5 (lower) for q = 1. Lines 
displayed by crosses are instability lines and line with black circles is Potts transition line. The 
model solid is thermodynamically unstable outside the instability lines. 
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averaging time of the order of 10^ MC steps per spin. 

One purpose here is to test the renormahzation group prediction of section III for the cases 
when the heat capacity of the Potts model is divergent at criticahty: q = 2,3,4. Though 
the renormahzation group analysis is exact for diamond hierarchical lattice, its predictions 
for the square lattice are to be checked. The other goal of the simulations is to learn about 
the influence of the discontinuous Potts transition for g > 4 on the equation of state of our 
model. To achieve this we simulate the g = 10 on the square lattice. 

To verify the accuracy of our simulations we compared the critical temperature estimated 
from our simulations to the exact critical temperature. For r = 1 it is: 

_ 1 + ^0 , . 

With q = 4, Eq = —0.5, we find Tc ^ 0.721. The Monte Carlo simulations give the same 
result up to 4 digits. 

In MC simulations, we work at finite sizes, so for each size we have to determine the 
"pseudo" transition which corresponds in general to the maximum of the specific heat or of 
the susceptibility. The maxima of these quantities need not to be at the same temperature. 
Only at the infinite size, they should coincide. The theory of finite-size scaling permits to 
deduce properties of a system at its thermodynamic limit. We have used in this work a size 
large enough to reproduce the bulk transition temperature up to the fourth decimal. We 
define the Potts order parameter Q by 

Q = [gmax(Qi, Q^, ..., Q,) - l]/{q - 1) (22) 

where Q„(ra = 1, ..., g) is the averaged value defined by 

Qn =< E '^(^J- - n)/{N,Ny) > (23) 

j 

Sj being the Potts spin at the site j. 

B. Results for q=10 

We recall that the Potts model shows a first-order transitions^ for g > 4 at a finite 
temperature. This is seen in Fig. |H] where the averaged Potts energy U defined by Eq. fllSp 
and the order parameter Q are shown for three values of r. 
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FIG. 8: Averaged Potts energy U and order parameter Q vs temperature T for q = W, r = 0.96 



(black circles), 1 (void circles), and 1.04 (black triangles), with N^ 
Lines are guides to the eye. 



N„ 



100 and E, 







-0.5. 



Note that these quantities show a large discontinuity at the transition temperature. These 
results confirm the first-order character of the transition. Repeating the simulations for other 
values of r, we determine the Potts transition line in the space {T,r) which agrees up to 
four digits with the exact critical line of Eq. (20), using Wc = vTO. Note that unlike in 
some first-order transitions, the slow heating and slow cooling of the system do not result 
in a hysteresis. The energy barrier between the two phases is believed to be therefore not 
so high. 

We calculate the stress a using the Eq. (1161) with u obtained from MC simulations shown 
above in Fig. [HI E{r) and dE{r)/dr being given by Eqs. ([T]) and (jH]). This is done for 
many values of r and T around the Potts transition curve in search for unstable regions 
predicted by the RG analysis shown in Figs. [THTl In practice, we fixed a temperature and 
then changed the value of r across the Potts transition line by following a vertical line in 
Fig. [71 In doing so, we obtained for each T the stress as a function of r. At a given T, if 
there is no crossing of the Potts transition line then the stress behaves as shown in the curve 
for T = 0.46 in Fig. [HI c goes smoothly through a minimum at a compression position to a 
maximum at a dilatation one. The solid is stable in the region between the minimum and 
the maximum since the modulus da/dr is positive. On the other hand, when the system 
crosses the transition line by varying r, the stress exhibits jumps as seen in the three curves 
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FIG. 9: Stress vs r for g = 10 at T = 0.36 (upper left), 0.38 (upper right), 0.40 (lower left), 0.46 



(lower right), with Nx 
comments. 
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-0.5. Lines are guides to the eye. See text for 



at T = 0.36, 0.38, 0.40 in Fig. |9l The discontinuity in stress is due to the Potts energy 
discontinuity associated with the first-order transition of the 10-state Potts model. The 
first-order nature in the case g = 10 observed above enhances the instabihty. As will be 
seen below, a qualitatively different instability occurs for g = 2,3,4 where the transition is 
of second order. 

It is interesting to show now in Fig. [TU]the stress versus T at a given r. As seen the stress 
undergoes a discontinuity at the Potts transition temperature if r 7^ 1, but the discontinuity 
vanishes for r = 1. 

To close this sub-section, we emphasize that for the system studied here, the size effects 
are indistinguishable from A^^. = A^ = 60 up. 

C. Results for q = 2, 3, 4 



Let us consider the case where g = 4. We show in Fig. [TT] the averaged Potts energy 
[Eq. flT5]) ] and the order parameter Q obtained by MC simulations as described above for 
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FIG. 10: Stress vs T for q = 10, at three values of r: cl.04 (black triangles), 1 (void circles) and 



0.96 (black circles), with N^ 
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-0.5. Note the discontinuities in the cases 



0.96 and 1.04. Lines are guides to the eye. See text for comments. 



q = 4, Eq = —0.5 and three values of r. The Potts transition is continuous with diverging 
heat capacity^". As a consequence of Eq. (10) the modulus becomes negative close to the 
transition temperature. This instabilty implies the emergence of a van der Waals loop in 
the dependence of the stress on r. We show this in Fig. [121 Hence the transition becomes 
weakly first order. 

For g = 3 and q = 2, the results are shown in Figs. [T3]fT4l and Figs. [T5lfT6l respectively. 
The replacement of the Potts continuous transition with a weakly first-order transition occurs 
here also because the Potts heat capacity is divergent for q = 2,3,4. As seen, the van der 
Waals loops, though weak, are present. 



V. CONCLUSIONS 



In this paper, we have used renormalization group and Monte Carlo simulations to study 
the phase diagram of a two-dimensional solid by using a model in which the interaction 
between neighboring atoms follows the Lennard- Jones potential. We have mapped the model 
into a g-state Potts model and investigated the effect of both a uniform compression and 
a uniform expansion of the volume of the solid. In the temperature, inter-atomic distance 
plane we find a line of Potts transitions and stability boundaries where the stress as a 
function of inter-atomic distance has an extremum. For the cases where the Potts heat 
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FIG. 11: Averaged Potts energy U and order parameter Q vs temperature T for q = 4, r = 0.96 
(black circles), 1 (void circles), 1.04 (black triangles), with Nx = Ny = 100 and Eq = —0.5. 

capacity is divergent (g = 2,3,4 for the square lattice, q > 6.8 for the hierarchical diamond 
lattice) a van der Waals loop (instability) occurs close to the Potts transition and thus a 
weak first-order transition (using Maxwell construction) replaces the continuous transition. 
This is a remarkable result that warrants further analysis. Monte Carlo simulations for 
large q values (where the square lattice Potts transition is discontinuous) indicate that the 
discontinuity in the Potts energy translates into a discontinuity in the stress as a function 
of the interatomic distance at the transition temperature. The influence of inter-atomic 
distance fluctuations on the nature of the phase transitions in this model will be addressed 
in a future study. Finally, let us note that while it is known that the solid phase cannot 
survive at finite T in two dimensions with continuous degrees of freedom of atom motions, 
our present study with discrete degrees of freedom shows some interesting behaviors which 
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FIG. 12: Stress (upper curve) and zooms (lower curves) vs r for g = 4, T = 0.67931, with 
N^ = Ny = 60 and Eq = -0.5. 



would serve as a starting point to study three-dimensional solids where melting mechanisms 
are not well understood.-*^ 

M.K. wishes to thank the University of Cergy-Pontoise for hospitality while this work 



20 



1.2- 



u 



1.4 



1.6- 



1.8 



1 • ■ ■ A * 



/0 






.:i^ 



J'^ 




FIG. 13: Averaged Potts energy U and order parameter Q vs temperature T for q = 3, r = 0.96 
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was carried out. 
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FIG. 16: Stress vs r near the instability region (upper) for q = 2 with several temperatures from 
T = 0.535 to 0.675: data points for T = 0.65 are marked with void squares. This case is shown for 
a large region of r in the lower figure. N^ = Ny = 60, Eq = —0.5. 



